Exercises and objectives

The objectives of the exercises of this assignment are:
1) Download and organise the data and model and plot staircase responses based on fits of logistic functions
2) Fit multilevel models for response times
3) Fit multilevel models for count data

REMEMBER: In your report, make sure to include code that can reproduce the answers requested in the exercises below (MAKE A KNITTED VERSION)
REMEMBER: This assignment will be part of your final portfolio

Exercise 1

Go to https://osf.io/ecxsj/files/ and download the files associated with Experiment 2 (there should be 29).
The data is associated with Experiment 2 of the article at the following DOI https://doi.org/10.1016/j.concog.2019.03.007

Put the data from all subjects into a single data frame

datadir<-"/Users/Laura/Desktop/Methods 3/github_methods_3_class/week_03/experiment_2"

laurafiles<-list.files(datadir,pattern='.csv',full.names=TRUE)

df <- vroom(laurafiles)

## Rows: 18131 Columns: 17

## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (5): trial.type, task, target.type, obj.resp, subject
## dbl (12): pas, trial, jitter.x, jitter.y, odd.digit, target.contrast, target...

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Describe the data and construct extra variables from the existing variables
1. add a variable to the data frame and call it correct (have it be a logical variable). Assign a 1 to each row where the subject indicated the correct answer and a 0 to each row where the subject indicated the incorrect answer (Hint: the variable obj.resp indicates whether the subject answered “even”, e or “odd”, o, and the variable target_type indicates what was actually presented.

obj.resp [i,16], o, e target_type [i,11], odd, even correct [i,19]

df$correct <- ifelse((df$obj.resp == "o" & df$target.type=="odd") | (df$obj.resp == "e" & df$target.type=="even"), 1,0)

describe what the following variables in the data frame contain, trial.type, pas, trial, target.contrast, cue, task, target_type, rt.subj, rt.obj, obj.resp, subject and correct. (That means you can ignore the rest of the variables in your description). For each of them, indicate and argue for what class they should be classified into, e.g. factor, numeric etc.

trial.type: a character variable encoded as either “staircase” (i.e. data was collected as a part of the adaptive staircase procedure which was performed at the beginning of the study) or “experiment” (i.e. data was collected as a part of the actual experiment trials). It should be re-encoded a factor in order for trial.type to be modelled later in the analysis, if relevant. However, if we were only interested in the experimental trials, we might consider filtering out all staircase trials*.

df$trial.type <- as.factor(df$trial.type)

pas: indicates one of the 4 categorically different ratings of the Perceptual Awareness Scale (PAS: Ramsøy & Overgaard, 2004), rated by the participant. It indicates the participant’s perceptual awareness of the stimulus where the numbers represent the following:*

(1) No experience (i.e. "no impression of the simulus. All answes are seen as mere guesses)*

(2) Weak Glimpse (i.e. "A feeling that something has been shown. Not characterized by any content, and this cannot be specified any further)*

(3) Almost Clear Experience (i.e. “Ambiguous experience of the stimulus. Some stimulus aspects are experienced more vividly than others. A feeling of almost being certain about one’s answer”)*

4) Clear experience (i.e. “Non-ambiguous experience of the stimulus. No doubt in one’s answer”)*

It should be encoded as a factor as it is a case of ordinal data.

df$pas <- as.factor(df$pas)

trial: refers to the trial number. It is numeric but should be factor

df$trial <- as.factor(df$trial)

target.contrast: indicates the contrast of the target stimulus, relative to the background, adjusted to match the threshold of each individual by using the QUEST-algorithm. It is numeric and should stay numeric as we are dealing with a continuous variable with values falling within the range of 0-1*

cue: cue shown to participants to indicate the set of possible digits that might appear on each trail. There were 36 different combination of cues (0-35). The task-column specifies the number of potential targets cued. As the cue variable contains nominal data, it should be reencoded into a factor.

df$cue <- as.factor(df$cue)

task: indicates the number of potential targets that were cued. Consisted of either 2,4, or 8 digits, indicated by the character variable task taking one of the 3 values:

singles (2 possible targets) e.g. 2:9 pairs (4 possible targets) e.g. 24:57 quadruplet (8 possible targets) e.g. 2468:3579

This variable should also be re-encoded as a factor as we are, again, dealing with nominal data.

df$task <- as.factor(df$task)

target_type: indicates the parity of the target stimulus (i.e. whether the digit was even or odd). It is encoded as a character but should be re-encoded as a factor.

df$target.type <- as.factor(df$target.type)

rt.subj: response time of answering Perceptual Awareness Scale Encoded as numeric as it is continuous data.

rt.obj: response time (time taken by participant to indicate the parity of the target stimulus). Encoded as numeric as it is continuous data.

obj.resp: the participant’s answer when asked to indicate the parity of the target stimulus (e = even, o = odd). Encoded as numeric but should be re-encoded as a factor.

df$obj.resp <- as.factor(df$obj.resp)

subject: participant ID, enabling us to distinguish between data collected from different participants. Is a character variable but should be re-encoded as a factor as we are dealing with a repeated measures design and therefore need to be able to model subject as a random effect in the analysis phase.

df$subject <- as.factor(df$subject)

correct: indicates accuracy of the respondent’s answer of the target stimulus’ parity. It is a binary variable where 0 indicates an incorrect response and 1 indicates a correct one, and should therefore be re-encoded as a factor.

df$correct <- as.factor(df$correct)

for the staircasing part only, create a plot for each subject where you plot the estimated function (on the target.contrast range from 0-1) based on the fitted values of a model (use glm) that models correct as dependent on target.contrast. These plots will be our no-pooling model. Comment on the fits - do we have enough data to plot the logistic functions?

library(tidyverse)
# make subset of df with only trial.type == staircase
df_staircase <- subset(df, df$trial.type == "staircase")

df_staircase$subject = gsub("(?<![0-9])0+", "", df_staircase$subject, perl = TRUE)
df_staircase$subject = as.integer(df_staircase$subject)
nopoolfun <- function(i){
  dat <- df_staircase[which(df_staircase$subject == i),]
  model <- glm(correct~target.contrast, family = 'binomial', data=dat)
  fitted <- model$fitted.values
  plot_dat <- data.frame(cbind(fitted,'target.contrast'=dat$target.contrast))
plot <- ggplot(plot_dat, aes(x = target.contrast, y = fitted))+
  geom_point()+
  geom_line(aes(x = target.contrast, y = fitted))+
  xlab('Target Contrast')+
  ylim(c(0,1))+
  theme_minimal()
print(plot)
}
# Running the function for each participant
for (i in 1:29){
  nopoolfun(i)
}

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

subjects <- c(1:16)
plots <- lapply(subjects, FUN=nopoolfun)

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

do.call(grid.arrange,  plots)

subjects <- c(17:29)
plots <- lapply(subjects, FUN=nopoolfun)

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

do.call(grid.arrange,  plots)

#  --------------- No pooling model  --------------- 
m_nopool <- glm(correct~target.contrast*subject, family=binomial(link=logit), data=df_staircase)
fitted_nopool <- fitted(m_nopool)
df_staircase$fitted_nopool <- fitted_nopool
# one plot 
ggplot(data = df_staircase)+
  geom_point(aes(x = target.contrast, y = fitted(m_nopool), color = subject))+
  geom_line(aes(target.contrast, fitted(m_nopool))) +
  labs(title = "No pooling plot")+
  facet_wrap(.~subject)

From visual inspection of the no pooling plot, it becomes clear that the fits are very bad. A subset data frame containing only the staircasing data doesn’t contain enough data to plot the logistic functions (the fitted values do not follow a sigmoid function).

on top of those plots, add the estimated functions (on the target.contrast range from 0-1) for each subject based on partial pooling model (use glmer from the package lme4) where unique intercepts and slopes for target.contrast are modelled for each subject

# Telling the model to expect differing baseline-levels of "correct" (the intercept, represented by 1) as well as differing responses to the main factor in question, which is "target.constrast" in this case.
m_partial <- glmer(correct~target.contrast + (1+target.contrast|subject), data = df_staircase, family = "binomial")
fitted_partial <- fitted(m_partial)
df_staircase$fitted_partial <- fitted_partial
## Plotting on top
ggplot(data = df_staircase)+
  geom_point(aes(x = target.contrast, y = fitted(m_nopool), color = "no pooling"))+
  geom_point(aes(x = target.contrast, y = fitted(m_partial), color= "partial pooling"))+
  geom_line(aes(target.contrast, fitted(m_partial), color = "partial pooling")) +
  geom_line(aes(target.contrast, fitted(m_nopool), color = "no pooling")) +
  labs(title = "No pooling and partial pooling model plot")+
  facet_wrap(.~subject)

in your own words, describe how the partial pooling model allows for a better fit for each subject

Partial pooling allows for a better fit for each subject as it takes into account the expected individual baseline differences depending on the stimulus contrast level as well as allows for different slopes for each subject. In other words, it allows the model to reflect that individuals would perform differently in different experimental settings. Partial pooling provides us a compromise between complete and no-pooling models, as it allows us to model both an average (general tendency of the whole dataset) and each level of the categorical predictor, subject. In partial pooling (a mixed effects model), we model both an average and each level, i.e. combine information from the population (fixed) effects of the complete pooling model and the subject-specific (random) effects of the no-pooling one. (((accounts for nested structure of our data, it assumes hierarchy in our data structure)))

Exercise 2

Now we only look at the experiment trials (trial.type)

Pick four subjects and plot their Quantile-Quantile (Q-Q) plots for the residuals of their objective response times (rt.obj) based on a model where only intercept is modelled

df_experiment <- df %>% filter(trial.type=="experiment")
df_experiment

## # A tibble: 12,528 x 18
##    trial.type pas   trial jitter.x jitter.y odd.digit target.contrast
##    <fct>      <fct> <fct>    <dbl>    <dbl>     <dbl>           <dbl>
##  1 experiment 1     0     -0.322     0.395          5          0.0392
##  2 experiment 1     1      0.499     0.495          3          0.0392
##  3 experiment 1     2      0.461    -0.198          5          0.0392
##  4 experiment 1     3      0.200    -0.0263         3          0.0392
##  5 experiment 1     4     -0.347    -0.0448         5          0.0392
##  6 experiment 1     5      0.0944   -0.106          5          0.0392
##  7 experiment 1     6      0.291    -0.347          5          0.0392
##  8 experiment 1     7      0.00541  -0.168          3          0.0392
##  9 experiment 2     8     -0.393    -0.357          5          0.0392
## 10 experiment 1     9     -0.110     0.185          3          0.0392
## # … with 12,518 more rows, and 11 more variables: target.frames <dbl>,
## #   cue <fct>, task <fct>, target.type <fct>, rt.subj <dbl>, rt.obj <dbl>,
## #   even.digit <dbl>, seed <dbl>, obj.resp <fct>, subject <fct>, correct <fct>

df_experiment$subject = gsub("(?<![0-9])0+", "", df_experiment$subject, perl = TRUE)
df_experiment$subject = as.integer(df_experiment$subject)
# pick 4 random from data frame
library(dplyr)
set.seed(1)
sample <- sample_n(df_experiment,4) # note that you could get the same number twice
sample # 03, 019, 012, 025

## # A tibble: 4 x 18
##   trial.type pas   trial jitter.x jitter.y odd.digit target.contrast
##   <fct>      <fct> <fct>    <dbl>    <dbl>     <dbl>           <dbl>
## 1 experiment 1     152      0.362   0.479          7          0.0628
## 2 experiment 1     227      0.450   0.188          5          0.0357
## 3 experiment 1     22       0.161  -0.269          9          0.106 
## 4 experiment 4     0        0.399   0.0117         5          0.0309
## # … with 11 more variables: target.frames <dbl>, cue <fct>, task <fct>,
## #   target.type <fct>, rt.subj <dbl>, rt.obj <dbl>, even.digit <dbl>,
## #   seed <dbl>, obj.resp <fct>, subject <int>, correct <fct>

# make df subsets for each subject
sub3 <- df_experiment %>% filter(subject==3)
sub19 <- df_experiment %>% filter(subject==19)
sub12 <- df_experiment %>% filter(subject==12)
sub25 <- df_experiment %>% filter(subject==25)
# make modes (lm cause continuous variables)
exp_model <- lm(rt.obj~1, data = df_experiment)   # overall
# FModels
m3 <- lm(rt.obj~1, data = sub3)
m19 <- lm(rt.obj~1, data = sub19)
m12 <- lm(rt.obj~1, data = sub12)
m25 <- lm(rt.obj~1, data = sub25)
# QQ-plot
qq3 <- ggplot(sub3, aes(sample = residuals(m3))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subject 003: Q-Q Plot for residuals") +
theme_bw()
qq19 <- ggplot(sub19, aes(sample = residuals(m19))) + stat_qq() +
stat_qq_line(colour = "red") + 
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subject 019: Q-Q Plot for residuals") +
theme_bw()
qq12 <- ggplot(sub12, aes(sample = residuals(m12))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subject 012: Q-Q Plot for residuals") +
theme_bw()
qq25 <- ggplot(sub25, aes(sample = residuals(m25))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subject 025: Q-Q Plot for residuals") +
theme_bw()
ggpubr::ggarrange(qq3, qq19, qq12, qq25)

## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4

i. comment on these
All 4 qq-plots indicates major deviances from normality, i.e. right skewness. We could try to improve the QQ-plots via log-transformation of the response time data.

does a log-transformation of the response time data improve the Q-Q-plots?

m3log <- lm(log(rt.obj)~1, data = sub3)
qq3log <- ggplot(sub3, aes(sample = residuals(m3log))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subj 3: QQ-plot residuals (log)") +
theme_bw()
m19log <- lm(log(rt.obj)~1, data = sub19)
qq19log <- ggplot(sub19, aes(sample = residuals(m19log))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subj 19: QQ-plot residuals") +
theme_bw()
m12log <- lm(log(rt.obj)~1, data = sub12)
qq12log <- ggplot(sub12, aes(sample = residuals(m12log))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subj 12:QQ-plot residuals") +
theme_bw()
m25log <- lm(log(rt.obj)~1, data = sub25)
qq25log <- ggplot(sub25, aes(sample = residuals(m25log))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subj 25: QQ-plot residuals") +
theme_bw()
ggpubr::ggarrange(qq3log, qq19log, qq12log, qq25log)

Doing a log-transformation improves the QQ-plots slightly, altough the data still deviates from normality.

Now do a partial pooling model modelling objective response times as dependent on task? (set REML=FALSE in your lmer-specification)

partial1 <- lmerTest::lmer(log(rt.obj)~task + (1|subject) + (1|trial), data = df_experiment, REML=FALSE)
partial2 <- lmerTest::lmer(log(rt.obj)~task + (1|trial), data = df_experiment, REML=FALSE)
partial3 <- lmerTest::lmer(log(rt.obj)~task + (1|subject), data = df_experiment, REML=FALSE)
partial4 <- lmerTest::lmer(log(rt.obj) ~ task + (1 + task | subject) + (1 | trial), data = df_experiment, REML = F)
partial5 <- lmerTest::lmer(log(rt.obj) ~ task + (1 + task | subject), data = df_experiment, REML=FALSE)

## boundary (singular) fit: see ?isSingular

## Warning: Model failed to converge with 1 negative eigenvalue: -7.6e+02

partial_model <- c("m1", "m2", "m3", "m4", "m5")
AIC <- c(AIC(partial1), AIC(partial2), AIC(partial3), AIC(partial4), AIC(partial5))
sigma <- c(sigma(partial1),sigma(partial2),sigma(partial3), sigma(partial4), sigma(partial5))
# calculating pseudo R-squared for mixed effects model (R2m = marginal R2, provides the variance explained only by fixed effects. R2c = conditional R2, provides the variance explained by the entire model, i.e., both fixed effects and random effects)
pacman::p_load(MuMIn)
r.squaredGLMM(partial1) # R2c = .2629

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.

##              R2m       R2c
## [1,] 0.006697664 0.2628567

r.squaredGLMM(partial2) # R2c = .0360

##              R2m        R2c
## [1,] 0.006663833 0.03597408

r.squaredGLMM(partial3) # R2c = .2253

##              R2m       R2c
## [1,] 0.006538142 0.2253338

r.squaredGLMM(partial4) # R2c = .2750

##              R2m     R2c
## [1,] 0.006742469 0.27504

r.squaredGLMM(partial5) # R2c = .2371

##              R2m       R2c
## [1,] 0.006538065 0.2370859

pseudo_r2 <- c(0.2629, 0.0360, 0.2253, 0.2750, 0.2371)
  
as_tibble(cbind(partial_model, sigma, AIC, pseudo_r2))

## # A tibble: 5 x 4
##   partial_model sigma             AIC              pseudo_r2
##   <chr>         <chr>             <chr>            <chr>    
## 1 m1            0.767353686788801 29459.9382110624 0.2629   
## 2 m2            0.877443929282292 32560.1503149819 0.036    
## 3 m3            0.786510836280231 29685.3126188651 0.2253   
## 4 m4            0.761030254763385 29327.6815811647 0.275    
## 5 m5            0.780526731035549 29560.3359954543 0.2371

which would you include among your random effects and why? (support your choices with relevant measures, taking into account variance explained and number of parameters going into the modelling)

m4 has the lowest sigma as well as AIC value. This is the model that predicts objective response times dependent on task (fixed effects), modelling a random intercepts for subjects and trial as well as random slopes for task. Random intercepts were modeled for subject as one wouold expect individual baseline performance differences among different participants. Trial as also been modelled as random intercepts as you would expect different trials to vary in difficulty (and participants might modify their performance as they get acquainted with the experiment).

explain in your own words what your chosen models says about response times between the different tasks

summary(partial4)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: log(rt.obj) ~ task + (1 + task | subject) + (1 | trial)
##    Data: df_experiment
## 
##      AIC      BIC   logLik deviance df.resid 
##  29327.7  29409.5 -14652.8  29305.7    12517 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.5363  -0.4958  -0.0279   0.5119   8.0799 
## 
## Random effects:
##  Groups   Name           Variance Std.Dev. Corr       
##  trial    (Intercept)    0.029952 0.17307             
##  subject  (Intercept)    0.149162 0.38621             
##           taskquadruplet 0.004001 0.06325   0.35      
##           tasksingles    0.032270 0.17964   0.38 -0.58
##  Residual                0.579167 0.76103             
## Number of obs: 12528, groups:  trial, 432; subject, 29
## 
## Fixed effects:
##                Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)    -0.26407    0.07317 29.75753  -3.609  0.00111 ** 
## taskquadruplet -0.07225    0.02054 28.60783  -3.518  0.00147 ** 
## tasksingles    -0.17868    0.03739 29.12028  -4.778 4.65e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tskqdr
## taskqudrplt  0.099       
## tasksingles  0.279 -0.114

The estimates for quadruplet and singles are negative, which signifies that these tasks have a shorter reaction time compared to the pair task. Relative to the quadruplet the singles task is generally faster.

The coefficients “taskquadruplet” and “tasksingles” are both negative and accompanied by small significant p-values, indicating that as the task changes from pairs to taskquadruplet and Tasksingles, respectively, the response time significantly decreases. The effect is, however, strongest from changing from pairs to singles (the respodents generally complete the singles tasks faster, relative to the quadruplet tasks).

Now add pas and its interaction with task to the fixed effects

# fails to converge: pastask <- lmer(log(rt.obj)~pas*task + (1 + task | subject) + (1 | trial), data = df_experiment, REML=FALSE) 
pastask <- lmer(log(rt.obj)~pas*task + (1 | subject) + (1 | trial), data = df_experiment, REML=FALSE)

how many types of group intercepts (random effects) can you add without ending up with convergence issues or singular fits?

pastask1 <- lmer(log(rt.obj) ~ task*pas + (1|subject), data = df_experiment, REML = FALSE)
pastask2 <- lmer(log(rt.obj) ~ task*pas + (1|trial) + (1|subject), data = df_experiment, REML = FALSE)
pastask3 <- lmer(log(rt.obj) ~ task*pas + (1|trial) + (1|subject) +  (1|target.contrast), data = df_experiment, REML = FALSE)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?

pastask4 <- lmer(log(rt.obj) ~ task*pas + (1|trial) + (1|subject) +  (1|cue), data = df_experiment, REML = FALSE)
pastask5 <- lmer(log(rt.obj) ~ task*pas + (1|trial) + (1|subject) +  (1|cue) + (1|target.contrast), data = df_experiment, REML = F)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?

The first 3 models (up to 3 group intercept) can be modelled without running into convergence issues if I add trial, subject and cue. However, if I try to add target.contrast as my third group intercept, I run into convergence issues. However, if adding cue as my third group intercept and target.contrast at my fourth, I can add 3 types (but run into issues at the fourth one).

create a model by adding random intercepts (without modelling slopes) that results in a singular fit - then use print(VarCorr(<your.model>), comp='Variance') to inspect the variance vector - explain why the fit is singular (Hint: read the first paragraph under details in the help for isSingular)

print(VarCorr(pastask5), comp='Variance')

##  Groups          Name        Variance 
##  trial           (Intercept) 0.0252761
##  cue             (Intercept) 0.0025684
##  target.contrast (Intercept) 0.0644419
##  subject         (Intercept) 0.0930045
##  Residual                    0.5592559

help(isSingular)

in your own words - how could you explain why your model would result in a singular fit?

We can see that the variances of all the random intercepts in the model pastask5 are all very low (especially our random-effect variance estimate for “cue” is nearly zero). Variances being estimated as 0 results in singular fit. It indicates that the model is overfitting, so the structure of our random effects is too complex for our data to support it.

Exercise 3

Initialise a new data frame, data.count. count should indicate the number of times they categorized their experience as pas 1-4 for each task. I.e. the data frame would have for subject 1: for task:singles, pas1 was used # times, pas2 was used # times, pas3 was used # times and pas4 was used # times. You would then do the same for task:pairs and task:quadruplet

pacman::p_load(dplyr)
data_count <- df %>% 
  group_by(subject, task, pas) %>% 
  dplyr::summarize("count" = n())

## `summarise()` has grouped output by 'subject', 'task'. You can override using the `.groups` argument.

data_count

## # A tibble: 340 x 4
## # Groups:   subject, task [87]
##    subject task       pas   count
##    <fct>   <fct>      <fct> <int>
##  1 001     pairs      1       109
##  2 001     pairs      2        45
##  3 001     pairs      3         4
##  4 001     pairs      4        12
##  5 001     quadruplet 1       141
##  6 001     quadruplet 2        23
##  7 001     quadruplet 3         3
##  8 001     quadruplet 4         1
##  9 001     singles    1        43
## 10 001     singles    2        87
## # … with 330 more rows

# the cause of this error is R’s confusion which summarize function (dplyr vs. plyr) it should use.

Now fit a multilevel model that models a unique “slope” for pas for each subject with the interaction between pas and task and their main effects being modelled

countmodel1 <- glmer(count~pas*task + (1+pas|subject), data = data_count, family = poisson, control = glmerControl(optimizer="bobyqa"))

summary(model1) i. which family should be used?

Poisson Regression should be used as it is best for modeling a dependent variable that consists of “count data” given one or more independent variable. Poisson distributed data is integer-valued, discrete, not continuous, and is limited to non-negative values. It also assumes that the errors follow a Poisson distribution rather than a normal distribution, and models the natural log of the response variable, ln(Y), as a linear function of the coefficients (rather than modeling Y as a linear function of the regression coefficients).

why is a slope for pas not really being modelled?

pas is encoded as a factor, meaning that it computes the analysis for each level seperately . Thus, we’ll get a slope for each separate level of pas, where each slope is relative to the reference level.

if you get a convergence error, try another algorithm (the default is the Nelder_Mead) - try (bobyqa) for which the dfoptim package is needed. In glmer, you can add the following for the control argument: glmerControl(optimizer="bobyqa") (if you are interested, also have a look at the function allFit)

The control argument “glmerControl(optimizer=”bobyqa“)” was added to avoid convergence error

when you have a converging fit - fit a model with only the main effects of pas and task. Compare this with the model that also includes the interaction

# Model with only the main effects of pas and task
countmodel2 <- glmer(count ~ pas + task + (1 + pas|subject), data = data_count, family = "poisson", control = glmerControl(optimizer="bobyqa")) 
rownames <- c("With interaction", "Without interaction")
AIC <- c(AIC(countmodel1), AIC(countmodel2))
resvar <- c(sum(residuals(countmodel1)^2),sum(residuals(countmodel2)^2)) #residual variance
res_sd<- c(sqrt((sum(residuals(countmodel1))^2/length(residuals(countmodel1))-2)),sqrt((sum(residuals(countmodel2))^2/length(residuals(countmodel2))-2))) # residual standard deviation 
as_tibble(cbind(rownames, AIC, resvar, res_sd))

## # A tibble: 2 x 4
##   rownames            AIC              resvar           res_sd           
##   <chr>               <chr>            <chr>            <chr>            
## 1 With interaction    3148.44094783359 699.486557245994 0.847368170996978
## 2 Without interaction 3398.54853775285 962.462811611566 1.32245901291909

The model that includes the interaction has the lowest AIC value as well as the lowest residual variance and residual standard deviation (as compared to the model fitted with only the main effects of pas and task).

indicate which of the two models, you would choose and why

I would choose the model that includes the interaction. Partly because of the lower AIC value, residual variance and residual standard deviation, but also because it is easier to theoretically justify choosing this model. We are not really interested in predicting count from PAS rating (i.e. the participants confidence rating). Rather, we are interested in the interaction between task and pas rating, as PAS is related to the specific task at hand and we would expect PAS rating to depend on the task, considering the fact that the task variable indicates the number of potential targets that were cued (consisted of either 2,4, or 8 digit) and thus also indicates the difficulty of the task (and a more difficult task would be expected to generate a lower confidence rating).

based on your chosen model - write a short report on what this says about the distribution of ratings as dependent on pas and task

summary(countmodel1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: count ~ pas * task + (1 + pas | subject)
##    Data: data_count
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   3148.4   3232.7  -1552.2   3104.4      318 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3871 -0.7853 -0.0469  0.7550  6.5438 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev. Corr             
##  subject (Intercept) 0.3324   0.5765                    
##          pas2        0.3803   0.6167   -0.75            
##          pas3        1.1960   1.0936   -0.84  0.63      
##          pas4        2.3736   1.5407   -0.86  0.42  0.72
## Number of obs: 340, groups:  subject, 29
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          4.03570    0.10976  36.769  < 2e-16 ***
## pas2                -0.02378    0.11963  -0.199 0.842458    
## pas3                -0.51365    0.20718  -2.479 0.013166 *  
## pas4                -0.77292    0.29075  -2.658 0.007853 ** 
## taskquadruplet       0.11490    0.03127   3.674 0.000239 ***
## tasksingles         -0.23095    0.03418  -6.756 1.42e-11 ***
## pas2:taskquadruplet -0.11375    0.04605  -2.470 0.013508 *  
## pas3:taskquadruplet -0.20901    0.05287  -3.954 7.70e-05 ***
## pas4:taskquadruplet -0.21500    0.05230  -4.111 3.94e-05 ***
## pas2:tasksingles     0.19536    0.04830   4.045 5.23e-05 ***
## pas3:tasksingles     0.24299    0.05369   4.526 6.02e-06 ***
## pas4:tasksingles     0.56346    0.05101  11.045  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pas2   pas3   pas4   tskqdr tsksng ps2:tskq ps3:tskq
## pas2        -0.742                                                     
## pas3        -0.829  0.613                                              
## pas4        -0.847  0.412  0.703                                       
## taskqudrplt -0.151  0.138  0.080  0.057                                
## tasksingles -0.138  0.126  0.073  0.052  0.484                         
## ps2:tskqdrp  0.102 -0.198 -0.054 -0.039 -0.679 -0.328                  
## ps3:tskqdrp  0.089 -0.082 -0.125 -0.034 -0.592 -0.286  0.402           
## ps4:tskqdrp  0.090 -0.083 -0.048 -0.093 -0.598 -0.289  0.406    0.354  
## ps2:tsksngl  0.098 -0.188 -0.052 -0.037 -0.342 -0.708  0.490    0.203  
## ps3:tsksngl  0.088 -0.080 -0.124 -0.033 -0.308 -0.637  0.209    0.486  
## ps4:tsksngl  0.092 -0.085 -0.049 -0.091 -0.324 -0.670  0.220    0.192  
##             ps4:tskq ps2:tsks ps3:tsks
## pas2                                  
## pas3                                  
## pas4                                  
## taskqudrplt                           
## tasksingles                           
## ps2:tskqdrp                           
## ps3:tskqdrp                           
## ps4:tskqdrp                           
## ps2:tsksngl  0.205                    
## ps3:tsksngl  0.184    0.451           
## ps4:tsksngl  0.507    0.474    0.427

data_count$task <- relevel(data_count$task, ref = "singles")
ggplot(data_count, aes(x = pas, y = count)) +
  geom_point(aes(pas, count), color = "blue") +
  facet_wrap(~ task) +
  theme_bw()

bp <- ggplot(data_count, aes(x=pas, y=count, group=pas)) + 
  geom_boxplot(aes(fill=pas))
bp

bp + facet_grid(task ~ .)

# bp + facet_grid(pas ~ task, margins=TRUE)

A general linear mixed effect model was fitted with a poisson link function to investigate the distribution of ratings (count) as dependent on pas and task. In other words, the model predicted our outcome variable count based on the fixed effects of pas and task as well as their interaction term. For random effects, pas was modelled as a random slope for each subject, with subject modeled as random intercept. The direction of the slope (whether it is positive or negative) is highly dependent on the interaction between pas and task.

When inspecting the plots above, it seems as if the distribution of ratings across the 3 tasks (i.e, pairs, quadruplet, singles) are quite similar (with PAS 3 being the least frequent rating across all taks). Nonetheless, it seems as if the singles task had a higher count of PAS ratings 4 (i.e., “Clear experience”) which one would expect as this is the least difficult task. From just eyeballing the boxplot, it becomes more difficult to notice the differences between pairs and quadruplets tasks.

include a plot that shows the estimated amount of ratings for four subjects of your choosing

# Using the same 4 respondents as previously (3, 19, 12, 25)
data.count_four <- data_count %>% 
  filter(subject == "003" | subject == "019" | subject == "012" | subject == "025")
m_four <- glmer(count ~ pas * task + (pas|subject), data = data.count_four, family = poisson)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0166871 (tol = 0.002, component 1)

data.count_four %>% 
    ggplot() +
      geom_point(aes(x = pas, y = fitted(m_four), color = "Estimated")) + 
      geom_point(aes(x = pas, y = count, color = "Observed")) +
      facet_wrap( ~ subject)

Finally, fit a multilevel model that models correct as dependent on task with a unique intercept for each subject

finalmodel<- glmer(correct ~ task + (1|subject), family = "binomial", data = df_experiment)
summary(finalmodel)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ task + (1 | subject)
##    Data: df_experiment
## 
##      AIC      BIC   logLik deviance df.resid 
##  13636.3  13666.0  -6814.1  13628.3    12524 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7450 -1.0599  0.4791  0.6483  0.9795 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  subject (Intercept) 0.3622   0.6018  
## Number of obs: 12528, groups:  subject, 29
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     1.11896    0.11774   9.503  < 2e-16 ***
## taskquadruplet -0.07496    0.05081  -1.475  0.14014    
## tasksingles     0.16603    0.05218   3.182  0.00146 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tskqdr
## taskqudrplt -0.220       
## tasksingles -0.213  0.494

does task explain performance?

Yes, task significantly explains performance (as measured by correctness) for all 3 task levels (p < .05 in all cases). In the quadruplets task, the performance is significantly worse as compared to the task pairs, but subjects’ performance in the singles task is better as compared to pairs task. This would be expected as an increase in task complexity would be expected to result in a decrease in the amount of correct answers.

add pas as a main effect on top of task - what are the consequences of that?

finalmodel2 <- glmer(correct ~ task + pas + (1 | subject), data = df_experiment, family = 'binomial')

Task looses is significance in explaining performance. Pas seems to be a better predictor (it significantly predicts correct). This makes sense as perceptual awareness theoretically must be believed to pressupose giving the correct answer (above the threshold that would be expected by chance alone.)

now fit a multilevel model that models correct as dependent on pas with a unique intercept for each subject

finalmodel3 <- glmer(correct ~ pas + (1 | subject), data = df_experiment, family = 'binomial')

finally, fit a model that models the interaction between task and pas and their main effects

finalmodel4 <- glmer(correct ~ task*pas + (1 | subject), data = df_experiment, family = 'binomial')

describe in your words which model is the best in explaining the variance in accuracy

rownames <- c("M1", "M2", "M3", "M4")
AIC <- c(AIC(finalmodel), AIC(finalmodel2), AIC(finalmodel3),AIC(finalmodel4))
resvar <- c(sum(residuals(finalmodel)^2),sum(residuals(finalmodel2)^2),sum(residuals(finalmodel3)^2),sum(residuals(finalmodel4)^2)) #residual variance
res_sd<- c(sqrt((sum(residuals(finalmodel))^2/length(residuals(finalmodel))-2)),sqrt((sum(residuals(finalmodel2))^2/length(residuals(finalmodel2))-2)),sqrt((sum(residuals(finalmodel3))^2/length(residuals(finalmodel3))-2)),sqrt((sum(residuals(finalmodel3))^2/length(residuals(finalmodel3))-2)),sqrt((sum(residuals(finalmodel4))^2/length(residuals(finalmodel4))-2))) # residual standard deviation 
as_tibble(cbind(rownames, AIC, resvar, res_sd))

## Warning in cbind(rownames, AIC, resvar, res_sd): number of rows of result is not
## a multiple of vector length (arg 1)

## # A tibble: 5 x 4
##   rownames AIC              resvar           res_sd          
##   <chr>    <chr>            <chr>            <chr>           
## 1 M1       13636.2881319358 13503.973068995  15.6567154698802
## 2 M2       12260.4130198953 12137.5692924277 13.6364434946894
## 3 M3       12256.8814750344 12138.0725449922 13.6379998755294
## 4 M4       12265.7719796063 12130.7523038036 13.6379998755294
## 5 M1       13636.2881319358 13503.973068995  13.6265684844832

anova(finalmodel, finalmodel2, finalmodel3,finalmodel4)

## Data: df_experiment
## Models:
## finalmodel: correct ~ task + (1 | subject)
## finalmodel3: correct ~ pas + (1 | subject)
## finalmodel2: correct ~ task + pas + (1 | subject)
## finalmodel4: correct ~ task * pas + (1 | subject)
##             npar   AIC   BIC  logLik deviance     Chisq Df Pr(>Chisq)    
## finalmodel     4 13636 13666 -6814.1    13628                            
## finalmodel3    5 12257 12294 -6123.4    12247 1381.4067  1     <2e-16 ***
## finalmodel2    7 12260 12312 -6123.2    12246    0.4685  2     0.7912    
## finalmodel4   13 12266 12362 -6119.9    12240    6.6410  6     0.3553    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

practical_exercise_3, Methods 3, 2021, autumn semester

Laura Emilie Hollænder Jensen

04/09-2021

Exercises and objectives

Exercise 1

Exercise 2

Exercise 3